% vim: set spelllang=en :

\documentclass[10pt]{article}

\usepackage[latin1]{inputenc}
\usepackage{amsmath, amssymb, amsfonts, amsthm} 
\usepackage{upgreek} 
\usepackage{amsthm} 
\usepackage{fullpage}
\usepackage{graphicx}
\usepackage{cancel}
\usepackage{wrapfig}
\usepackage{subfigure}
\usepackage{mathrsfs}
\usepackage{outlines}
\usepackage[font={sf,it}, labelfont={sf,bf}, labelsep=space, belowskip=5pt]{caption}
\usepackage{hyperref}
\usepackage{titling}
\usepackage{xifthen}
\usepackage{color}
\usepackage{fancyhdr}
\usepackage[title]{appendix}
\usepackage{float}
\usepackage{enumitem}

\usepackage{bm}
\usepackage{minted}
\usepackage{xcolor}

\pagestyle{fancy}
\headheight 24pt
\headsep    12pt
\lhead{\documenttitle}
\rhead{\today}
\fancyfoot[C]{}
\lfoot{}
\rfoot{\thepage}
\renewcommand{\headrulewidth}{0.4pt}
\renewcommand\footrulewidth{0.4pt}
\newcommand{\documenttitle}{Assignment Solution}
\newcommand{\coursetitle}{SE2324: Mathematical Foundation of Computer Sciences(Spring 2021)}
\newcommand{\authorname}{Student Name, ID1001001001} % Your name and ID here

\makeatletter

\setlength{\droptitle}{-50pt}

\title{\documenttitle \vspace{1ex} \\ \Large \coursetitle \vspace{-1ex}}
\author{\authorname\vspace{-1ex}}


% BEGIN DOCUMENT
\begin{document}
\maketitle


%%%%%%%%%%%%%%%%%%%
\section{README}

... 

%%%%%%%%%%%%%%%%%%%
\section{NumPy Warm-up (15 points)}

Consider the following one-dimensional \textit{Gaussian probability distribution}, also known as the Normal distribution or bell curve distribution:

\begin{center}
$
G(x \mid \mu, \sigma)=\frac{1}{Z} \exp \left[-\frac{1}{2} \frac{1}{\sigma^{2}}(x-\mu)^{2}\right] \quad
$
where 
$
Z=\sqrt{\frac{\pi}{\frac{1}{2} \frac{1}{\sigma^{2}}}}
$
\end{center}

Here, the function $G : \mathbb{R} \rightarrow \mathbb{R}$ takes as input a scalar $x$, and produces as output a scalar. 
The particular bell shape of $G$ is determined by two parameters: the mean $\mu \in \mathbb{R}$ and the variance $\sigma^2 \in \mathbb{R}$

Numerically verify the following identity in Python:

\begin{center}
$\int_{\mathbb{R}} G(x) d x=1$
\end{center}

\begin{enumerate}[label=2.\arabic*]
    \item (10 points)
    Choose a few different one-dimensional Gaussian functions (by choosing different mean and variance values), plot them.
    \item (10 points)
    Verify the above identity for each Gaussian function. 
\end{enumerate}



%%%%%%%%%%%%%%%%%%%
\section{Numerics and Linear Algebra (60 points)}\label{sec:num_lin_alg}

In this question, you will explore how the condition number of a matrix can have practical influence on which algorithms (e.g. LU, Cholesky) you can use to solve linear systems of the form: $A\vec{x}=\vec{b}$ where $A\in \mathbb{R}^{n\times n}$.
We will study the Vandermonde matrix, as well as a matrix constructed from a finite Fourier Series basis.
In this question, you will interpolate samples of the analytical function

\begin{center}
$f(x)=\frac{1}{1+x^2}$    
\end{center}

for $x\in[0,1]$. Your monomial interpolants (with $N=n+1$ terms) are given by

\begin{center}
$g_V(x)=\sum_{j=0}^{n} c_j x^j$
\end{center}

and

\begin{center}
$g_{F}(x)=\sum_{j=1}^{N / 2} c_{j} \sin (j \pi x)+\sum_{j=N / 2+1}^{N} c_{j} \cos ((j-N / 2) \pi x)$.
\end{center}

Use $M=m+1$ uniformly sampled positions

\begin{center}
$x_i=ih,\text{ }i=0\dots m$,
\end{center}

with spacing $h=1/m$, to generate $M$ samples of the test function $f$,

\begin{center}
$f_i=f(x_i),\text{ }i=0\dots m$.
\end{center}

To estimate the polynomial coefficients $\vec{c}$, you will assemble and solve the linear systems

\begin{center}
$V\vec{c}=\vec{f}$
\end{center}

and

\begin{center}
$F\vec{c}=\vec{f}$
\end{center}

where $V$ is the M-by-N Vandermonde matrix with entries

\begin{center}
$V_{ij}=(x_i)^j$
\end{center}

and $F$ is the finite Fourier Series basis matrix with entries

\begin{center}
$F_{i, j-1}=\left\{
    \begin{array}{ll}
        \sin \left(j \pi x_{i}\right), & \text { if } 1 \leq j \leq N / 2 \\
        \cos \left((j-N / 2) \pi x_{i}\right), & N / 2+1 \leq j \leq N
    \end{array}
\right.$
\end{center}

For simplicity, assume $M=N$ for the entire problem.

\begin{enumerate}[label=3.\arabic*]
    \item (25 points) \label{q:3.1}
    Use an LU solve (\textit{scipy.linalg.lu} from SciPy package) to estimate the monomial coefficients $\vec{c}$. 
    Report the residual L2 norm for both linear systems when $N = 8$ and $N = 16$. 
    
    Hint: After LU decomposition, you can solve the resulting linear systems step by step using the NumPy API \textit{numpy.linalg.solve}().
    
    Hint: Given a problem $A\vec{x}=\vec{b}$, the residual L2 norm is $||r||_2=||A\vec{x}-\vec{b}||_2$.
    
    \item (10 points)
    Using the \textit{numpy.linalg.cond} function in NumPy, plot $N$ vs. $cond(V)$ and $N$ vs. $cond(F)$ for $N = 4, 6, 8, ...32$. 
    Write a couple of sentences explaining the reasons for the trends in these two plots. 
    
    Hints: Use a logarithmic scale in y axis for better clarity. Also, try wrapping the creation of $V$ and $F$ into functions that you can call repeatedly to generate the required output data. These functions will be helpful for the next part.
    
    \item (15 points) \label{q:3.3}
    A necessary condition for being able to use Cholesky factorization is that the matrix must be positive definite. 
    Construct $A_V = V^TV$ and $A_F = F^TF$ for $N = 4, 6, \dots32$ (a total of 30 matrices). 
    Mathematically, when would these matrices be positive definite? Explain.
    Using the \textit{isposdef}() function introduced in Appendix \ref{app:posdef}, check to which matrices NumPy reports as positive definite. Create a table of values that includes the following columns: $N$, \textit{isposdef}($A_V$), \textit{isposdef}($A_F$), cond($V$), cond($F$).
    What is the largest value of $N$ where $A_V$ is positive definite, and what is the condition number of that $V$? 
    What is the largest value of $N$ where $A_F$ is positive definite, and what is the condition number of that $F$? 
    Are these condition numbers connected in some way? If so, how?
    
    \item (10 points)
    For $N = 8$, transform the linear systems above into symmetric positive definite systems according to Question \ref{q:3.3}, and use Cholesky factorization (\textit{numpy.linalg.cholesky}()) to solve them. Report the residual L2 norm for each solution. Compare the residuals to Question \ref{q:3.1}: how does Cholesky compare to LU?


\end{enumerate}


%%%%%%%%%%%%%%%%%%%
\section{Least Squares Problems and QR (25 points)}\label{sec:lin_lsp}

In the question above, a $M\times M$ square linear system is solved for the interpolation coeffcients.
The resulting interpolation function can provide a solution that pass through all the sample points.
However, in some applications, finding such a solution could be too time-consuming and unnecessary, or even impossible (Consider two sample with the same $x$ value and different $y$ values).
In such cases, we usually reduce the number of the interpolants (let $M>N$), and solve the resulting over-determined linear system using least square method.
Instead of looking for the exact solution to an over-determined system, we look for the solution with the smallest error vector $V\vec{c}-\vec{f}$ (or $F\vec{c}-\vec{f}$) by minimizing the overall backward error:

\begin{center}
minimize $\|V \vec{c}-\vec{f}\|_{2}^{2}$
\end{center}

\begin{enumerate}[label=4.\arabic*]
    \item (15 points)
    Solve the least square system with QR decomposition(\textit{numpy.linalg.qr}()) when $M=16$, $N=4,8$.
    \item (10 points)
    Plot the $g_V$, $g_F$ when $M=16$, $N=4,8$, compare them with the analytical function $f(x)$ and the interpolation function obtained in Question \ref{q:3.1}.
\end{enumerate}


\end{document}
